######################################################################
#Replication file for "The Impact of China's AIIB on the World Bank"
#Authors: Jing Qian, James Vreeland, Jianzhi Zhao
#Codes to replicate Figure A.8
######################################################################
rm(list = ls())

#To install the bpCausal package (ver 0.0.1)
# install.packages('devtools', repos = 'http://cran.us.r-project.org') # if not already installed
# devtools::install_github('liulch/bpCausal')

#Load bpCausal package
library(bpCausal)

#Other packages used
library(tidyverse)
library(RColorBrewer)

#Setwd
setwd("C:/Users/qianj/Dropbox (Princeton)/AIIB_WB_Replication")

#Load custom functions
#Note: Given the number of tests performed in the appendix, 
#several custom functions are created to implement estimation methods in a more efficient way.
#These functions merely create a wrap-up of existing functions, like bpCausal and gsynth,
#in order to avoid typing parameters repeatedly.

source("code/custom_functions.R")

###################################
#Estimation for Figure A.8
###################################
#----------------
#(1) Load data
#----------------
load("data/data_main.RData")
load("result/bpcausal_original.RData") #Original results

#--------------
#(2) Setup
#--------------
dv.list = c("hard_project_count_all")
covs.all = c("gdppc_log_lag",
             "population_log_lag",
             "election_either_lag",
             "fdi_gdp_lag",
             "debt_gni_lag",
             "oda_gni_lag",
             "polity2_lag",
             "unsc_lag",
             "IdealPointDistance_lag")
X.use = Z.use = A.use = covs.all

#---------------------------------------
#(3) Estimation: Exclude non-founding members
#~3.3 minutes
#---------------------------------------
#Exclude non-founding members joined by 2019
iso2c.nonfounder.2019 = df.main %>%
  subset(aiib_member_year <= 2019) %>%
  subset(aiib_non_founder == 1) %>%
  .$iso2c %>%
  unique()

df.use = df.main %>%
  subset(!iso2c %in% iso2c.nonfounder.2019)

D.use = "aiib_founder_2016"

result.a8.1 = bpcausal.group(dv.list = dv.list,
                             df.use = df.use,
                             D.use = D.use,
                             X.use = X.use,
                             Z.use = Z.use,
                             A.use = A.use)

#Number of observations
df.use %>%
  .[, c("iso2c", "year", dv.list, D.use, covs.all)] %>%
  .[complete.cases(.),] %>%
  nrow()
  

#---------------------------------------
#(4) Estimation: Treatment = Non-founding members
#~2.5 minutes
#---------------------------------------
#Exclude founding members
df.use = df.main %>%
  subset(aiib_founder == 0)

D.use = "aiib_non_founder_vary"

result.a8.2 = bpcausal.group(dv.list = dv.list,
                             df.use = df.use,
                             D.use = D.use,
                             X.use = X.use,
                             Z.use = Z.use,
                             A.use = A.use) 

#Number of observations
df.use %>%
  .[, c("iso2c", "year", dv.list, D.use, covs.all)] %>%
  .[complete.cases(.),] %>%
  nrow()

###################################
#Produce Figure A.8
###################################
#-------------
#(1) Prepare results
#-------------
#Summarize
result.use = effSummary(result.original)
result.founder.only.use = effSummary(result.a8.1$hard_project_count_all)
result.nonfounder.only.use = effSummary(result.a8.2$hard_project_count_all)

#-----------------
#(2) Parameters
#-----------------
#Parameters
n = 1
xlim = c(-2, 2)
cex = 2
xtext = 0.6
colors = brewer.pal(n = 3,
                    name = "Dark2")
gap = 1

#-------------------
#(3) Plot
#-------------------
#----------
#Setup
png(filename = "figure/Figure_A8.png",
    height = 800,
    width = 1200)

par(mar = c(5, 1, 5, 1))

#----------
#Empty plot
plot(1,
     type = "n",
     xlab = "Estimated Average Treatment Effect on the Treated",
     ylab = "",
     yaxt = "n",
     xaxt = "n",
     xlim = xlim,
     ylim = c(0, 4),
     cex.lab = 2)
box()

#------------------
#Hard
#------------------
#Original
points(x = result.use$est.avg[1],
       y = 3,
       col = colors[1],
       pch = 19,
       cex = cex)

arrows(x0 = result.use$est.avg[2],
       x1 = result.use$est.avg[3],
       y0 = 3,
       y1 = 3,
       length = 0.05,
       angle = 90,
       code = 3,
       col = colors[1],
       cex = cex,
       lwd = cex)

#Excluding non-founding members
points(x = result.founder.only.use$est.avg[1],
       y = 2,
       col = colors[2],
       pch = 19,
       cex = cex)

arrows(x0 = result.founder.only.use$est.avg[2],
       x1 = result.founder.only.use$est.avg[3],
       y0 = 2,
       y1 = 2,
       length = 0.05,
       angle = 90,
       code = 3,
       col = colors[2],
       cex = cex,
       lwd = cex)

#Exclude founding members
points(x = result.nonfounder.only.use$est.avg[1],
       y = 1,
       col = colors[3],
       pch = 19,
       cex = cex)

arrows(x0 = result.nonfounder.only.use$est.avg[2],
       x1 = result.nonfounder.only.use$est.avg[3],
       y0 = 1,
       y1 = 1,
       length = 0.05,
       angle = 90,
       code = 3,
       col = colors[3],
       cex = cex,
       lwd = cex)
#------------------
#Other elements
#------------------
#Abline
abline(v = 0,
       lty = "dashed",
       lwd = cex,
       col = "gray")
#Legend
legend("topright",
       legend = c("Original Result",
                  "Exclude Non-founding members",
                  "Treatment = Non-founding members"),
       cex = cex,
       col = colors[1:3],
       lty = c(1, 1),
       bty = "n",
       lwd = cex)

#Axis
axis(side = 1,
     at = seq(min(xlim),
              max(xlim),
              1),
     cex.axis = cex)

#Title
title("World Bank Infrastructure Projects",
      cex.main = cex,
      line = 1)
#----------
#Close
dev.off()